Anomalous scaling in passive scalar advection: Monte-Carlo Lagrangian trajectories 
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We present a numerical method which is used to calculate anomalous scaling exponents of structure 
functions in the Kraichnan passive scalar advection model (R. H. Kraichnan, Phys. Fluids 11, 
945 (1968)). This Monte-Carlo method, which is applicable in any space dimension, is based on 
the Lagrangian path interpretation of passive scalar dynamics, and uses the recently discovered 
equivalence between scaling exponents of structure functions and relaxation rates in the stochastic 
shape dynamics of groups of Lagrangian particles. We calculate third and fourth order anomalous 
exponents for several dimensions, comparing with the predictions of perturbative calculations in 
large dimensions. We find that Kraichnan's closure appears to give results in close agreement with 
the numerics. The third order exponents are compatible with our own previous nonperturbative 
calculations. 



The Kraichnan rapid advection model Q is a simplified 
model for turbulent advection of a passive scalar T{r) 
in which the velocity field scales in space but decorre- 
lates in time infinitely rapidly. The scaling properties 
of this model are interesting: The structure functions 
S n (R) = ([T(r + R) - T(r)] n ) depend on R like power 
laws, S n (R) ~ R> n , and the exponents Cn are "multiscal- 
ing" in the sense that they are nonlinear functions of the 
index n in contradiction with classical (Kolmogorov) ap- 
proaches. The importance of this model lies in the fact 
that it has been the first example in which the mechanism 
for multiscaling, which appears also in Navier-Stokes tur- 
bulence ||], has been identified and in which a con- 
trolled calculation of scaling exponents is possible. As 
such it is an important laboratory for testing ideas that 
may pertain also for the understanding of multiscaling 
in Navier-Stokes turbulence. It is not surprising there- 
for that it drew enormous attention and theoretical ef- 
forts. To date, the scaling exponents in the Kraichnan 
model have been calculated only in a small part of the 
parameter space, which consists of the velocity scaling 
exponent < £ < 2 (see below), and the space dimen- 
sionality d. Perturbative calculations were performed in 
the nearly integrable limits £<lj§,2-£<l [|L|I, 
and d ^> 1 PUl2]]. In addition, Kraichnan proposed g] 
an interesting closure method to estimate the exponents 
for even n and arbitrary d. Direct numerical simulations 
were performed only in two dimensions [^],[l0| . The expo- 
nents £2 and £3 were computed throughout the parameter 
space, with an explicit solution in the former case, and 
by a finite difference scheme in the latter jll| . 

In a recent article |jL2f , we presented a new interpreta- 
tion of the phenomenon of anomalous scaling in passive 
scalar correlation functions: the theory connects between 
the scaling exponents £ n and the relaxation rates of the 
geometric shape formed by n points advected by the tur- 
bulent flow. Starting from a generic distribution in the 
space of shapes, the stochastic advection can be consid- 
ered as a relaxation to the invariant measure in this space, 



while the over-all scale undergoes a Richardson diffusion. 
In |Q this identification was used to perform l/d per- 
turbative calculations of the scaling exponents. In this 
Letter we apply the theory, using Monte-Carlo methods, 
to calculate the scaling exponent throughout the entire 
parameter space |13[ ] . We generate random sample paths 
for the n trajectories, and measure the shape relaxation 
rate by the standard method of analyzing the autocorre- 
lation of the signal. Compared to direct numerical simu- 
lation, this method demands relatively modest comput- 
ing resources: there is no need to keep track of the whole 
scalar field. The realizations of random velocity are gen- 
erated only at the instantaneous positions of the advected 
points, and this amounts to sampling a small set of cor- 
related Gaussian random variables. Hence, it is possible 
to perform simulations in arbitrary integer dimensions. 
Another advantage of the Monte-Carlo method is that 
it is trivially parallelizable, in the sense that the data is 
collected from numerous independent realizations, whose 
calculation is completely independent, and may there- 
fore be performed on different machines. Nevertheless, 
the accurate calculation of the anomalous exponents is 
still not easy, due to slow convergence rate and the large 
number of realizations that is therefore required. 

The Kraichnan model describes the advection of a pas- 
sive scalar by a random velocity field u which is Gaus- 
sian, 5-correlated in time, incompressible and self similar 
in space: 

((u(r,t)-u(r',*))® (u(r,t')-u(r',i'))> 
= h(r - r')d(t - 1') , 

The passive scalar is driven by a large scale forcing (of 
scale L), which is taken Gaussian and (5-correlated in time 
as well. The fundamental objects of the Kraichnan model 
are the many point correlation functions 

T 2n {v 1 ,v 2 ,...,v 2n ) = {T{r 1 ,t)T{r 2l t)...T{r 2n ,t)) , (2) 
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where pointed brackets denote an ensemble average with 
respect to the stationary statistics of the forcing and the 
statistics of the velocity field. 

The correlation functions may be expanded in a series 
of scale invariant terms, i.e., terms which are homoge- 
neous functions of their variables. The most important 
term for R/L — > is expected to be a zero mode of 
the Kraichnan partial differential operator operator |||| . 
The functions Ti n contain, in addition to their own zero- 
modes, also contributions from the zero modes of lower 
order functions and normal scaling contributions due to 
the external forcing. These do not contribute however 
to the nth order structure functions, and thus do not 
affect the scaling exponent ( n . The question is how to 
compute the homogeneity exponents of the zero modes. 
Kraichnan himself avoided this (difficult) problem and 
suggested a closure procedure whose main result is a 
closed form ansatz for the scaling exponents 



C2„ - \ 



(3) 



where (2 = 2 — £. We will compare the £ = 1 line of this 
prediction to our numerics. 

We explain now shortly how multiscaling of scalar cor- 
relation functions can be obtained from the random evo- 
lution of shapes as a function of the overall scale. More 
details can be found in p^ |. We consider n points ad- 
vected simultaneously by the random velocity field u, 



9tTi(t) =Jl(Ti(t),t), 



1 < i < n. 



(4) 



The configuration of these n points can be described 
at each instant by one overall scale s, and a set of di- 
mensionless variables Z which describe its geometry. For 
example, in the case n = 3, the overall scale could be de- 
fined by the sum of edge lengths, and the shape by two 
of the angles of the triangle defined by the three points. 
As far as scaling properties are considered the precise 
definition of the overall scale and the shape variables is 
unimportant, and see fl2]| for details. 

The evolution of the shape defined by our n advected 
particles is characterized by an operator that we denote 
by 7„(Zo, Z, -^). This operator determines the probabil- 
ity that, starting at scale Sq with shape Z , the n points 
will form the shape Z, for the first time that the overall 
scale reaches the value s > sq. We are interested in the 
asymptotics of the shape evolution operator; its large s 
expansion is 



7n(Zo, Z, • 



so 
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/3 m (Z) Mro (Z ). (5) 



/3o(Z) is the invariant measure of the shape dynamics 
(A n ,o = 0), and the functions /3 m (Z), m > are the 'ex- 
cited states' of these dynamics, describing the decay into 
the invariant measure, with increasing exponents A„. m . 



Ref. p2] established a correspondence which implies that 
the exponents A n , m , are precisely the set of anomalous 
scaling exponents of the nth order correlation function. 
The functions /x m 's are the adjoint family of left eigen- 
functions of 7„, and are the corresponding scaling struc- 
tures, namely the zero modes of the Kraichnan operator. 

The basic ingredient for the numerical implementation 
of the shape dynamics for the calculation of scaling expo- 
nents is a finite difference approximation of the solution 
of eqs. (^) describing the random Lagrangian trajectories 
of the n points. We used the simple Ito-Euler method for 
this purpose [OJ , which makes the approximation 

r i (t + *t)~r i (t)+Vft»h(t), (6) 

where the rji are a set of Gaussian random vectors whose 
covariance as a function of the rVs is given by the spatial 
part of the covariance of the n random vectors u(r^). 

The results of the difference scheme are then used 
to calculate the instantaneous scale s(t) using some con- 
venient definition of s (we used the square root of the sum 
of squares of point separations). We fixed a set of thresh- 
old values for s spaced logarithmically, and for each event 
where the value of s increases for the first time beyond 
one of the threshold values, we record the instantaneous 
value of a certain function of the shape c(Z). Once more, 
the function a is largely arbitrary; however, it has to be 
a permutation-symmetric function of the coordinates. 

To extract the relaxation rates, it is convenient to cal- 
culate the autocorrelation of the signal <x(s). In fact (for 
s < s'), 

{a(s)a(s')) = J dZdZ'(3 (Z) ln (Z,Z',s'/s)a(Z)a(Z') = 



(a) 2 + {<thx)o ( CT >i 



+ •••, (7) 



where () denote average with respect to f3 m . Thus, the 
autocorrelation function of a relaxes to a constant, and 
the asymptotic relaxation rate is equal to the nth order 
structure function scaling exponent £„. 

The precise measurement of the scaling exponent is not 
easy, however, since the signal contains a constant whose 
fluctuation tend to mask the decaying part which con- 
tains the scaling information. Additionally, in the case 
of n = 4 which we consider below, three point contri- 
butions contaminate the signal with relatively persistent 
transients. For this reason a very large number of realiza- 
tions, O(10 6 ), is required to measure a scaling exponent, 
and the function a needs to be tuned carefully to decrease 
the effect of transients p!E[ . Indeed, the autocorrelation 
for any generic definition of a decays asymptotically with 
an exponent £ n , however, for a poorly chosen definition 
of cr, the transients may dominate for a very large range. 
We have tuned our definitions for a by trial and error to 
expose the asymptotic decay rate as early as possible. 
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We removed the constant from the autocorrelation sig- 
nal using two methods, the first was by numerically differ- 
entiating the auto-correlation signal, and the second by 
subtracting the value of the autocorrelation at the far- 
thest separation measured (s/sq = 100). Both methods 
give consistent results though the numerical differentia- 
tion is more noisy while subtraction tends to introduce a 
bias in the large s results. 

In our n — 3 calculations we used 

0-3 = Iog(lia) + log(Zi 3 ) + log(Z 23 ) (8) 
and in the four point 

<T4 = (h^f + (W 24 ) 4 + (M23) 4 . (9) 



Where lij is the distance between points i,j normalized 
such that J2i<j l ij = !■ 

The case of third order correlation function exponents 
deserves attention, being the simplest non-trivial one. 
The Lagrangian paths are relatively easy to generate in 
this case, and the number of realizations needed to ob- 
serve scaling behavior is smaller than in the four point 
case. Additionally, demonstrated that three point 
correlation function can be calculated directly by a finite 
difference scheme, and thus the two calculations can be 
tested against each other for consistency. 

We generated a set of data for three point relaxation 
rates in 4 dimensions for various values of £. The num- 
ber of required realizations was O(10 5 ), except for small 
values of £ where a larger number is required. We also 
generated data for £ = 1 for several space dimensionali- 
ties, as in the four point case discussed below. In fig. 1 
we display a typical data set measuring the relaxation of 
the three point distribution to its invariant measure. 
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Similar results were obtained for the three point simu- 
lations for other parameter values, with generally a longer 
scaling range the closer is £ to 2. The results are summa- 
rized in fig. 2, along with the leading scaling exponents 
calculated by the method of O], We find satisfactory 
agreement between the predictions of the two methods, 
providing support for the validity of the present theoret- 
ical and numerical framework. 
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FIG. 2. The leading three point scaling exponent in four 
dimensions, calculated as a function of £ by the present 
Monte-Carlo method (circles) and by the method of [ll]is 
difficult to estimate the systematic errors in the results of the 
latter method, so no error bars are given. However, we have 
indications that these errors are smaller the statistical errors, 
except when £ is near 2. In the latter case we use results 
whose accuracy is higher than those published in [11] 

Possibly the main interest in our method is in its ability 
to provide measurement of £4 from the simulation of four 
point relaxation. In measuring £4 we concentrated on a 
scan through various dimensions and fixed £ = 1. The 
possibility of performing measurements in various dimen- 
sions (including the 'physical' d = 3 case) exemplifies the 
superiority of this algorithm with respect to direct nu- 
merical simulations, where dimensions larger than 2 are 
inaccessible due to memory requirements. The (i-scan al- 
lows us to test, for the first time, the perturbative large 
d expansion results [§],[l2). In fig. 3 we show the results 
of a typical four point calculation. 
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FIG. 1. The increments of (o"3(s')(J3(s / s)) as a function of 
s, for a three point simulation in 4 dimensions with £ = 1, 
normalized by the normal scaling value s 2 . Taking increments 
eliminates the constant contribution (0-3) 2 . The dashed line is 
a (shifted) least squares power law fit. The data are averaged 
over 2 x 10 5 realizations 
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FIG. 3. The increments of {(7i(s')(7i(s' s)) from a four point 
simulation, with the same parameters as in fig .1, normalized 
in the same way, with a power law fit (dashed line). The 
upward tilt of the graph means that the anomaly is positive. 
The data are averaged over 1.6 x 10 6 realizations 

The results of the three and four point simulations as 
a function of d are summarized in fig. 4, along with the 
prediction of large d perturbation theory, and the pre- 
diction of Kraichnan's closure. We find that Kraichnan's 
prediction agrees well with the (4 results, whereas the 
perturbative calculations are still not relevant. We can- 
not estimate how large d should be to agree with the 
large d perturbation theory. The prediction for £ 3 (given 
in [12]) is in rough agreement with our simulations. 
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FIG. 4. The leading three point (circles) and four point 
(squares) scaling exponent for velocity scaling exponent £ = I 
as a function of the inverse space dimensionality. The upper 
and lower dashed lines show the prediction of first order per- 
turbation theory in 1/d for the three and four point cases, 
respectively. The continuous line is Kraichnan's closure pre- 
diction Eq.(4) 

In summary, Lagrangian trajectory approach (devel- 



oped independently by p3| ) is seen to provide a numeri- 
cal tool allowing us to compute the anomalous exponents 
of the Kraichnan model, in principle for the entire param- 
eter space. The simulations are based on a new interpre- 
tation of anomalous scaling linking it to relaxation rates 
in the stochastic shape dynamics; the success of the simu- 
lations provides further support to this physically appeal- 
ing picture of anomalous scaling in this model. Finally 
we note that the basic principles employed here are not 
specific to the Kraichnan model, and it might be possible 
to use them in other contexts, perhaps sacrificing some 
of the simplicity. In particular, the fact that we need 
to consider only n points simultaneously depends heav- 
ily on the (5-correlated nature of the velocity, and would 
probably not carry over to more realistic situations. 
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We need approximately 20 cpu-days to obtain 10 6 itera- 
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tions on a typical alpha workstation. Using a cluster of 
8 such stations a point in parameter space can be calcu- 
lated relatively accurately within 2.5 days. 
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